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■ The energy splitting E^ a in two and four dimensional Ising models is measured in 

a cylindrical geometry on finite lattices. By comparing to exact results in the two 
dimensional Ising model we demonstrate that E^ a can be extracted very reliably from 
Monte Carlo calculations in practice. In four dimensions we compare the measured 
Eo a with two different theoretical predictions on the finite size behavior of the energy 
splitting. We find that our numerical data are in favor of the predictions based on the 
semiclassical dilute instanton gas approximation. 
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1 Introduction 



In studies of first order phase transitions on finite lattices one faces the phenomenon of 
tunneling between different states. The tunneling rate is related to the surface tension, 
which is the free energy per unit area associated with the interface separating these different 
states. As the surface tension is the fundamental parameter for the nucleation rate in first 
order phase transitions the determination of it has attracted a lot of attention. Examples are 
widespread. In d=2, models like the g-state potts model are of interest for statistical physics. 
In d=3 binary liquids give direct experimental information about the surface tension that can 
be compared to theoretical calculations using the Ising model as an effective description. In 
d=4 it is the first order QCD finite temperature phase transition where attempts have been 
made to extract the surface tension. This list of examples is of course by far not exhausting 
and many more examples could be found. 

Here we study the two and four dimensional Ising models below the critical temperature 
as examples of a typical system undergoing a first order phase transition between two ordered 
states. In a cylindrical geometry where the time dimension is much longer than the linear 
spatial dimension, the Ising system will break into domains along the time direction through 
tunneling. The inverse of the correlation length along the time direction defines a energy 
splitting between the ground state and the first excited state. This energy splitting is related 
to the surface tension U which, in turn, is associated with the domain wall energy. The 
finite volume dependence of the energy splitting has been a long standing problem. A 
few years ago, a definite finite size behavior for the energy splitting was derived based on 
a semiclassical dilute instanton gas approximation 0. The result agreed with an earlier 
intuitive conjecture by Brezin and Zinn- Justin ||. Furthermore, a relation between the 
surface tension and the renormalized parameters was derived in this approximation. This 
relation was checked in a numerical simulation and good agreement had been found both in 
four H and three dimensions ||. 

Recently, Borgs and Imbrie || gave a new derivation for the finite size behavior of the 
energy splitting for a general Ising like system. Their result is claimed to be exact. However, 
it is in disagreement with the instanton calculation. 

In this paper we confront both theoretical predictions with a Monte Carlo simulation. 
We first calculate the energy splitting in the two dimensional Ising model. By comparing 
to the known exact formula, we demonstrate that the method developed in ref. [§]] is very 
efficient and reliable for calculating the energy splitting. Then we turn our attentions to the 
four dimensional Ising model. 
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2 The Ising model 



The Ising model in d-dimensions is denned by 



where the fields $ can only assume values ±1. Here A is the volume of the lattice, i.e. TL 
in d=2 and TL 3 in d=4 with T the time and L the spatial extent of the lattice. Periodic 
boundary conditions are taken in all directions. The hopping parameter k is related to the 
inverse temperature. As is well known, the model eq. ([!]) can be solved exactly in d=2 
and has therefore played the role as a very often used toy model to test ideas and methods. 
In d=4 the Ising model is in the same universality class as the $ 4 theory. It exhibits therefore 
properties like triviality and was used to study at least qualitatively important aspects of 
the Standard Model physics §, g g|. 

The Ising model eq. (|]) is symmetric under parity transformation $ — > — $. For k > k c 
with k c = |log(l + V2) « 0.22036 in d=2 and k c « 0.0748 in d=4, this global Z(2) 
symmetry is spontaneously broken in infinite volume. We have a double well potential with 
minima at +M and —M. There are two degenerate ground states |0+ > and |0_ > which 
can be intuitively thought of as wavepackets centered around +M and — M, respectively. 
|0+ > and |0_ > have no overlap and a complete Hilbert space for the system can be built 
on either one of these ground states. Assuming we choose |0 + > as the ground state, the 
system will have a nonzero magnetization < + |$|0 + >= M. If one includes an external 
source j the magnetization will align with j. With varying external source the system will 
pass a first oder phase transition at j = and the magnetization M jumps from +M to 
—M. 

In a finite system the parity symmetry is never spontaneously broken. For k > k c the 
|0+ > and |0_ > states are no longer the ground states and have finite overlap. A wavepacket 
centered around one minimum of the effective potential can always tunnel and spread over 
to the other minimum. The ground state |0 S > is symmetric under parity transformation 
and may be expressed as 

|0 s >=-^=(|0 + >+|0_>). (2) 

The antisymmetric linear combination of |0+ > and |0_ > forms the parity odd first excited 
state |0 o > 

|0 a >=-^=(|0 + >-|0_>). (3) 

If we set the ground state energy to be zero, there is a small energy splitting E 0a between 
the first excited state |0 a > and the ground state \0 S >. This energy splitting is due to the 
tunneling between |0 + > and |0_ > states. 



As is well known JT|, [10], 11] the energy splitting E Qa vanishes exponentially fast with the 
spatial extent of the lattice and is related to the surface tension a, E Qa ~ exp(— crL^ 1 ). 
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The finite size behavior of Eo a can be analyzed by choosing a lattice geometry where the 
time extent T is much larger than the spatial extent L of the lattice. In this so called 
cylindrical geometry the system breaks up into domains of different magnetization, separated 
by interfaces in the spatial direction. 

Defining the timeslice operator St 

St= E ^> (4) 

for large enough t the tunneling energy can be obtained from the correlation function 

< S Q S t >= Cl (e- Eoat + e-^-*)) . (5) 



We calculated the correlation function < SoS t > numerically in the d=2 and d=4 Ising 
models in the broken phase. We used the cluster algorithm |12L |13[ in its 1-cluster version 
where we performed 10 updates in between measurements. We used both, the correlation 



function < SqSi > and its improved estimator version [13 



A 



<s„s t >=< n ^e(^iP)> . 



(6) 



x,y 



Here |c| is the size of the cluster C and is the cluster incidence function 



1 if both x,y G C 
else. 



(7) 



To achieve a realistic error for the energies E 0a we split the runs into several subruns and 
fitted the data according to eq. (|) for each subrun. We then determined the values of E Qa 
by the mean and its error of these individual results. 



3 The two dimensional Ising model 



We have chosen the two dimensional Ising model to test the method of extracting the surface 
tension from the finite size behavior of E 0a . In two dimensions the L-dependence of E 0a in 
the cylindrical geometry is known exactly |TU| , 



Ex 
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'2k' 



where 



7(x) = cosh 1 [cosh(2/t*) cosh(4/t) — cos(7nc)] 
and tt* is the dual coupling 

tanh(/t*) = 



-4k 



(8) 

(9) 
(10) 



3 



This exact expression has an asymptotic expansion for L —>■ oo 



2 sinh(cr) 



(i + o( z )) 



(11) 



Oa — 



7rL 



where the surface tension a is given by 



1 + e 
1 -e 



-4k "1 



a = 2(2k - «*) = 4k - In 



(12) 



Our lattices had a geometry LT where T = 1000 and L varied from 5 — 16. We have 
chosen two k values in the broken phase of the model, k = 0.2326 and k = 0.2252. We 
performed 50,000 measurements with 10 cluster updates in between each measurement. We 
fitted the spin-spin correlation function eq. (||) and the improved estimator eq. @ according 
to eq. (P) and performed a block error analysis as described above. Both operators gave 
consistent results. We checked that our fit results are stable against varying the start and 
end point of the fits. 

The resulting values of E$ a together with the exact formula eq. <M) is displayed in Fig. 1. 
We find a perfect agreement between the numerical data and the exact formula. This shows 
that the energy splitting E^ a can be extracted from the correlation function very precisely 
and reliably. 

We remark here that although E^ a can be measured accurately, extracting the surface 
tension from E^ a is not easy in this case. The surface tension is only clearly defined when 
L is large enough and E 0a can be well approximated by its asymptotic form eq. (]TT|). In 
d = 2 eq. ([TT|) sets in very slowly. For example, at k = 0.2326, E 0a calculated from eq. ( |i~l"D 
is 20% smaller than the exact value even when L = 16. In higher dimensions we expect the 
situation to improve since the asymptotic behavior of E 0a should set in much earlier with 
an error of 0(1/ L d ~ l ). 



4 The four dimensional Ising model 

It is well known that the Ising model is in the same universality class as the A$ 4 theory. 
Using the language of the continuous A$ 4 field the energy splitting E 0a can be calculated in 
a dilute instanton gas approximation [pl| . If we write the bare potential as 



v ; 4 24 ' 



(13) 



the energy splitting can be written as 




(14) 
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where det' is the determinant without the zero mode and the classical action S c comes from 
a single instanton configuration 

S c = 2^L d ~ 1 , (15) 
and the determinants M and Mo are given by 

M = -d 2 + m 2 - ^mlcosh' 2 (^4) , (16) 

and 

M = -d 2 + m 2 . (17) 

In eq. ([[4]) various factors can be understood intuitively. The factor exp(— S c ) comes from 
summing over a noninteracting multi-instanton configuration, the \f~S~ c factor comes from an 
integration over the location of the instantons and the ratio of determinants det'M/detM 
comes from integrating over the Gaussian fluctuations around the classical instanton config- 
uration. The classical instanton action eq. flT5] ) scales with the spatial volume L d ~ x . The 
ultraviolet part in the determinant ratio \det 'M j ' detM§\~ x l 2 will renormalize the bare pa- 
rameters in S c and the infrared part was conjectured to contribute an additional factor L^ 1 
0]. Thus the energy splitting has a general form 

Eoa = CinstL^e aL , (18) 

where a is the surface tension and u> is conjectured to be (d — 3)/2 under the dilute instanton 
gas approximation. 

Minister has calculated the determinant ratio det'M/detM in a finite spatial volume 
0]. The result of his work is in agreement with the conjectured form eq. ([Tj^) with u = 
(d— 3)/2. In addition the surface tension and the constant Ci nst were calculated in terms of 
the renormalized parameters. For example, in d = 4, we have 



3 

m. 



r 

9r 



9r I , 7T 



+ 0(g 2 ) 



(19) 




16tt 2 \8 4^3, 

and 

C inst = 1.650584/2^ . (20) 

V 9r 

In eqs. ( p!9|) and ( p0| ) m r and g r are the renormalized mass and renormalized quartic coupling 
respectively. These quantities are input parameters which have to be determined. A high 
statistic Monte Carlo simulation was performed at k — 0.076 [ffl. Eo a , m>r and g r were 
measured independently and the relation given by eq. (|19D was confirmed. 

Recently, Borgs and Imbrie claimed 0] to have calculated the exact finite size behavior 
of the energy splitting. In d — 4 they get 

E 0a = 2e- aLS , k > n c . (21) 
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This formula differs from eq. flI8| ) in two aspects. First, the prefactor in front of exp(— o~L 3 ) 
is independent of the volume. Second, the numerical value of the pre-exponential factor 
equals two and is believed to be independent of k. 

In order to confront both theoretical predictions with numerical simulation results we 
performed Monte Carlo simulations to determine E 0a . In the same manner as in d=2 we 
studied the correlation function eq. (|5]) and its improved estimator eq. (||) also in the four 
dimensional case. Here we have chosen four k values where the results for k = 0.076 are 
taken from reference |4[]. The measured E^ a values are listed in Table 1. As mentioned 
above the values of g r and m T are needed as input for the coefficient Ci ns t and the surface 
tension Oi nst in the instanton calculation. For the k values (beside k = 0.0760) used here we 
measured m r and A r on large lattices in a similar way as in [|J. The resulting values of m r 
and g r can be found in Table 2. 

We show the finite size behavior of E^ a for all four k values in Fig. 2. The solid lines are 
the predictions from the instanton calculation eqs. ([18]), (|T9|) and (|20|) and the solid circles 
are our measured E 0a . One of the assumptions in the instanton calculation is that the linear 
size of the spatial volume L is large. Thus we expect the agreement with MC data will only 
set in when L > 2£ with £ = l/m r being the correlation length in a hypercubic geometry. 
This seems to be the shown in Fig. 2. Our largest spatial volume is L = 14 at 

k = 0.0751 which is not obviously large enough so that the asymptotic behavior for E 0a has 
set in. At k — 0.0754 and L = 12 we checked that there is 18% finite size effect when T 
is increased from 120 to 160 and the effect is reduced to 6% at L = 11. This means that 
T should be at least more than 10 times larger than L to suppress the finite size effects in 
time direction. Thus an increase in L would require a large increase in T which will require 
a substantial increase in computer time. Given our limited numerical data, it is not possible 
to determine the power uj of L dependence in the pre-exponential factor in eq. (|I8"D . Instead, 
at each k value we take 3 data points with the largest L and fit to the form 

\n^ = -a flt L 3 + \nC flt . (22) 

The fitted results for af it and Cf it are listed in Table 3, along with x 2 per degrees of freedom. 
For comparison we also list in Table 3 the predicted a inst and C inst values from eqs. (0) and 
(p0|). It is clear from Table 3 that the measured Ofn and Cfu values are in reasonable agree- 
ment with the instanton calculation prediction. Given the fact that the one-loop correction 
term in eq. ([H]) contributes from 9% at k = 0.0751 to 13% at n = 0.0764, it is unlikely that 
the observed agreement with the fitted numerical results is accidental. 

We also tried a fit to the form 

E 0a = Ce- aLZ (23) 

and the results are listed in Table 4. From \ 2 both fits eq. fl2"2] ) and eq. ( |2"3"| ) are of the 
same quality with a slightly smaller \ 2 f° r ec L- (p3)- The fitted values for a are not very far 



from those listed in Table 3. However, the fitted pre-exponential factor C values are much 
smaller than 2 and are dependent on the value of k. This is a contradiction to eq. (pT|). 



A direct comparison to Borgs and Imbrie result is not possible because in eq. ( pT|) 
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the surface tension a is not expressed in terms of the renormalized parameters and should 
be treated as a free parameter. However, since both fits to eq. (p2|) and eq. (|23|) give 
approximately the same values for the surface tension as shown in Tables 3 and 4 and those 
values are in reasonable agreement with the predicted value a inst given by eq. (|19"D, we may 
take ai ns t value as an estimate of the true surface tension. If we use the value of ai nst in eq. 
( pT|) and setting the coefficient to 2 as predicted by Borgs and Imbrie, the calculated £"00 
will have a huge difference from MC data as shown in Fig. 3 (n = 0.0754). 



5 Conclusion 



Choosing a cylindrical geometry where the time extent of the lattice is much larger than 
the spatial extent, we studied the finite size behavior of the energy splitting E 0a on periodic 
lattices. We demonstrated in the two dimensional Ising model that the values of E 0a can be 
extracted from the spin-spin correlation function very precisely. The finite size dependence 
of E 0a was compared to the known exact results. We find perfect agreement between the 
numerical simulation results and the exact solution for E 0a as shown in Fig. 1. 

Extending this method to the four dimensional Ising model we confronted our numerical 
results with two theoretical predictions, one coming from a dilute instanton gas approxima- 
tion 0|, eq.(|18D and one recent work, claiming to have obtained the exact finite size behavior 
for E 0a ||, eq. (|21~1) . We measured the infinite volume values of m r and g r (as they are re- 
quired as input for the instanton calculation) at four different k values (see Table 2). The 
measured values of E Qa (see Table 1) agree reasonably with the predictions by eqs. (|T8l) , ( |i~9"D 
and (pCf ) from the dilute instanton gas approximation as shown in Fig. 2. Our fits to the 
form of eq. (|22|) or eq. ( p3|) give comparable \ 2 P er degrees of freedom. This suggests that 
our numerical data are not accurate enough to determine the power of the L dependence 
in the pre-exponential factor of E Qa . However, assuming E 0a has the form of eq. ([18]) with 
uj = (d — 3)/2, the fitted values for a and C agree well with the predictions given by eqs. 
( |19"D and (|20"D as shown in Table 3. In contrast, a fit to the form eq. (0), shown in Table 4, 
gives C values that are much smaller than two and k dependent. This is in contradiction to 
the result of Borgs and Imbrie, eq. fl2~T|). 



We conclude that given the choice between the instanton calculation and the result of 
ref. 0, our numerical data favor the instanton calculation result. Two comments are in 
order. First our numerical simulations are performed close to k c so that the system is inside 
the scaling region. The instanton calculation is only valid inside the scaling region since 
the domain walls are treated as continuous surfaces and renormalization language is used. 
Second, at present we can not determine the value of u in eq. ([Tj|) from our numerical 
data. Hypothetically, it is possible to have a scenario that E 0a takes the form in eq. ( p3|) 
with a k dependent C. Close to the critical point, C would be much smaller than two and 
gradually increases to two when k becomes large. We emphasize that so far this scenario has 
no theoretical backing. In order to determine the value of u further numerical simulation is 
needed on a larger lattice and with better statistics. 
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Table Caption 



Table 1: We give the values of the energy splitting Eo a as determined from fits to the 
correlation function eq. (^). The lattice geometry has been TL 3 where T = 120 for k = 
0.0760 and T = 160 elsewhere. The values at k = 0.0760 (except the L=ll result) are from 
reference [|J. 

Table 2: We give the renormalized mass m r and the renormalized coupling g r for the 
four dimensional Ising model for all k values used. The values for k = 0.0760 are taken from 
reference [|J. 

Table 3: We take E 0a measured on the largest three volumes at each k value and fit the 
data according to eq. (p2|) . The fitted afu and Cfu are listed here along with \ 2 per degrees 
of freedom. For comparison a inst and Ci nst calculated from eqs. (|19|) and ( p0|) are also listed. 

Table 4: Similar to Table 3, but data are fitted to the form of eq. (P3"|). Note C is much 
smaller than two and dependent on k. 



Figure Caption 

Figure 1: The energy splitting E^ a in the two dimensional Ising model as obtained from 
the Monte Carlo simulation. The open symbols are at k = 0.2326, the solid symbols are 
at k = 0.2252. The errorbars are smaller than the size of symbols. The lines represent the 
exact solution eq. @. Note that no fit is involved. 

Figure 2: The energy splitting E 0a for the four dimensional Ising model at four k values. 
The solid circles indicate the measured E 0a . The statistical errors are smaller than the size 
of symbols. The full curves are the predictions from the instanton calculation eq. ([18]) 0. 

Figure 3: The energy splitting E^ a for the four dimensional Ising model at k = 0.0754. The 
dashed curve is given by eq. (|21|) with a calculated from eq. (|i~9|). For comparison the full 
curve is from the instanton calculation eq. (fl8|) . 
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Table 1 
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Table 4 
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